Constrained Grid-Based Filter

ABSTRACT

A constrained four dimensional grid-based filter (CGBF) used to provide state estimates for a target maneuvering in two dimensions. Optimal grid and sampling sizes or chosen and the kinematic constraints of the target a y used to restrict the possible predicted states resulting in quality target estimates.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims the benefit of priority from U.S. ProvisionalPatent Application Ser. No. 61/708,673, filed on Oct. 2, 2012.

STATEMENT OF GOVERNMENT INTEREST

The invention described herein may be manufactured and used by or forthe Government of the United States of America for governmental purposeswithout payment of any royalties thereon or therefore.

BACKGROUND

The present invention relates to an improved system for using a gridbased filter to track a maneuvering target. More specifically, butwithout limitation, the present invention is a constrainedmulti-dimensional grid based filter used to provide state estimates fora maneuvering target.

Over the years, many systems have been developed to address the targettracking problem. However, tracking multiple maneuvering targets, asshown, in FIG. 1, remains a notable weakness of these systems. At theheart of a tracking system is the target track filter. The track filtercomputes the most likely target state given all the measurementscollected on that target. This most likely state is typically based onthe one that minimizes the mean squared error (MSE). In the early 1960s,Rudolf Kalman published a paper entitled “New Approach to LinearFiltering and Prediction Problems.” (Transactions of the ASME—Journal ofBasic Engineering, 82 (Series D): 35-45.) Kalman described a filter thatyields minimum MSE state estimates for systems that exhibit lineardynamics. This filter became known as the Kalman filter and since itsdebut, it has been extensively studied and has emerged as the premiermethod for target track filtering due to its computational efficiencyand its filtering effectiveness.

The Kalman filter is now used in many applications from GPS tracking toair traffic control. If the motion is linear and the measurement errorsare Gaussian, a Kalman filter can provide accurate estimates of a targetstate efficiently.

Nevertheless, although the Kalman filter is efficient and effective forcomputing state estimates of a moving target, it produces poor estimateswhen tracking a maneuvering target. The reason for these poor estimatesis that the Kalman filter is a linear filter. It assumes the dynamics ofthe tracked target can be modeled as a linear system, but maneuveringtargets do not generally exhibit linear motion. The reason the targetmotion needs to be linear is due to a property of the normaldistribution. A normal random variable remains normal through a lineartransform. If the transform is not linear, there is no guarantee thatthe normal variable will remain normal.

The Kalman filter includes a process noise parameter which provides ameans to capture the uncertainty in the target motion. Thus, astraightforward way to allow the Kalman filter to track maneuveringtargets would be to simply employ a large process noise. The idea is totreat the non-linear target motion simply as a larger uncertainty in itslinear motion. However, increasing the process noise parameter causeslarge errors in the computed state estimates. These larger errors causemore ambiguity making the tracking process more difficult.

Since the Kalman filter was so effective for linear targets, researchersbegan looking for enhancements and/or extensions to the Kalman filter.One of the early extensions was the extended Kalman filter (EKF). Radarwas (and still is) the dominant sensor used for target tracking. One ofthe early refinements to radar systems was the ability to measure theclosing rate of the target. This closing rate is known as the rangerate. Actually, the range rate is a closing rate only when the target ismoving radially towards, or away from the radar. If the target is movingalong some other path, the range rate is measuring only that componentof the target's velocity that is moving radially towards (or away from)the radar. This radial component is proportional to the cosine of thedifference of the angles between the target's course and the radialangle from the target to the radar. This radial angle to the target isoften called the azimuth angle. Cosine is a non-linear function andsince the range rate is a function of cosine, the (linear) Kalman filterhad to be extended to address this non-linear quantity. The resultingExtended Kalman Filter (EKF), includes this range rate measurement intothe calculations by performing a first-order approximation of itsadditional information. Since the Kalman filter is inherently a linearfilter, the intent is to linearize the range rate information. However,the cosine cannot always be estimated well with only a first-orderapproximation, so the EKF is a suboptimal filter. As a result, an EKFwill at times perform no better than a Kalman filter and could evenperform worse by failing to properly track the target.

Several other extensions to the Kalman filter subsequently emergedincluding iterated extended Kalman filters (IEKF), unscented Kalmanfilters (UKFs), and interacting motion model Kalman filters (IMMs). Allthese methods have the advantage that they are based on the Kalmanfilter and have closed-form solutions. This makes them computationallyefficient and well-grounded in a mathematical formulization. However,all these methods make assumptions about the maneuvering target that areoften not true, which leads to erroneous estimates.

Vast improvements to radar, computers and computational power have ledresearchers to pursue more numerical (and thus, more computational)methods for target tracking. Two methods that emerged were grid-basedfilters (GBFs) and particle filters (PFs). These methods were broadlyclassified as sequential Monte Carlo (SMC) techniques. A key feature ofthese methods is that they allow many of the restricted assumptionsplaced on the closed-form methods to be relaxed. It then became possibleto consider more general, non-linear target dynamics. Any presumedmotion of the target could be modeled without being limited to lineardynamics. In addition, the measurement errors could be modeled morerealistically, unlike the Kalman filter-like methods, which assumeGaussian error on the measurements even though radar measurements aregenerally not Gaussian (at least not in rectilinear space).

The idea behind a GBF is to first discretize the region where thetargets are presumed to reside into a set of cells. The next step is tomodel how the targets would transition from cell to cell over time.Finally, measurements are used to increase the probability of the targetbeing in a cell. Unfortunately, even with faster computers, thisapproach is far too computational to execute in real-time. Although thegrid-based filter approach is an attractive solution, due to itscomputational burden, it is difficult to implement and, thus, has notgained widespread use. What is needed is a computationally economicalmethod capable of using a four dimensional grid based filter fortracking a target maneuvering in two dimensions.

SUMMARY

The present invention is directed to a multi-dimensional grid basedfilter for tracking a moving target. The filter provides state estimatesfor a maneuvering target using a grid-based filter by initially forminga first grid of cells to capture an initial target state comprising anestimate of the target's position, course and speed using a firstmeasurement. A second grid of cells, of a finite size and scale, isformed wherein each cell comprises possible estimated states to whichthe target could have transitioned over a time between the first and asecond measurement. Using a sampling technique, the predictiondistribution is estimated for each cell in the grid using the variousmaneuvering probabilities of the target. Each cell corresponds to aparticular rectangular region in an x and y coordinate plane. An initialestimate is made of the probability that the target is in thecorresponding rectangular region. There is also an estimate of theexpected initial speed of the target and the variance of the speed inthe corresponding rectangular region. Similarly, there is an estimate ofthe expected initial course of the target and the variance of the coursein the corresponding rectangular region. A discretized probabilitydistribution over all the cells in the second grid is formed using MonteCarlo sampling. The position distribution is assumed to be uniformlydistributed over the cell and each cell forms a normal distribution forthe possible speeds for that region. Also, each cell forms a normaldistribution of the possible courses for that region. This processrepeats for each measurement update where the second grid becomes thefirst grid for the next update.

DRAWINGS

These and other features, aspects and advantages of the presentinvention will become better understood with reference to the followingdescription and appended claims, and accompanying drawings wherein:

FIG. 1 is a picture of an aircraft tracking various surface targets viaradar.

FIG. 2 is an algorithm for computing the destination grid (i.e., thesecond grid) for each update.

FIG. 3 is a graph of true predicted positional containment region.

FIG. 4 is a table comparing the Constrained Grid-Based Filter predictedstate and covariance to the Kalman Filter as the maximum turn a of thetarget is varied.

FIG. 5 is an illustration showing the grid cell update using Monte Carlosampling.

FIG. 6 is a graph shoving the target motion symmetry, mirrored statecalculation.

DESCRIPTION

In the following description of the present invention, reference will bemade to various embodiments which are not meant to be all inclusive. Thecurrent invention can be implemented using various forms of software forexecution on a variety of computer systems. The preferred embodiments ofthe present invention are illustrated by way of example below and in thereferenced Figures.

The power of a grid based filter (GBF) can be exploited withtwo-dimensional grid. The solution assumes that the grid discretizes thex,y locations, and then extends what is stored in each cell of the grid.Instead of merely storing the probability of the target being at a givenx,y location, four other values are also included: 1) the expected speedof the target given it is at that location, 2) the variance of speeds ofthe target at that location 3) the expected course of the target at thegiven location and 4) the variance of the targets expected course at thegiven location.

Because the kinematic constraints on a target are usually specified interms of speeds and turn rates, these values are used instead ofcomponent velocities. This ensures that the targets kinematicconstraints are enforced and exploited in the analysis. It is notnecessary to store the velocity as speed and course. It could be storedas x and y velocity components. The purpose of storing the course andspeed is to more easily and directly exploit the speed, acceleration,and turn rate constraints of the target. However, it should be notedthat working with course angle brings up different computationdifficulties due to its circular nature.

There are two grids, an originating grid and a destination grid. Thesegrids can be thought of as the From and To grids. Let the From grid belabeled F and the To grid be labeled T. Then F(x,y) is the cell atlocation x,y in the From grid. ‘Dot notation’ will be used to selecteach of the five values within a cell where F(x,y).prob, F(x,y).spd,F(x,y).svar, F(x,y).crs, and F(x,y).cvar are the probability, meanspeed, speed variance, mean course, and course variance, respectively.The algorithm in FIG. 2 shows how the probability, mean speed and speedvariance values are computed. The computation of the mean course and itsvariance require additional steps to deal with the circular nature ofthe course parameter, but are not further described.

There are three separate loops in FIG. 2. The first loop forms theweighted sums (probability) based on Monte Carlo samples. It is also theloop that determines to which track a new measurement should be assignedin multi-target scenarios. The measurement would be assigned to the gridwhich yields the highest cumulative probability. After the sampling iscomplete, the second loop converts these weighted sums into mean andvariance, using the equations 1-3 shown below. The final loopre-normalizes the mass in the cells.

$\begin{matrix}{P = {\sum p_{i}}} & {{Eq}.\mspace{14mu} 1} \\{\overset{\_}{s} = {\frac{1}{P}{\sum{p_{i}s_{i}}}}} & {{Eq}.\mspace{14mu} 2} \\{\sigma_{s}^{2} = {{\frac{1}{P}{\sum{p_{i}s_{i}^{2}}}} - \left( \overset{\_}{s} \right)^{2}}} & {{Eq}.\mspace{14mu} 3}\end{matrix}$

The 2D Target Motion Model

The algorithm in FIG. 2 also references three external functions:Normal, RandomMove, and CellOverlapFunction. The ‘Normal’ functioncreates normally distributed random variables. The RandomMove functioncreates random, but kinematically feasible predicted states over thetime interval. The function allows any target motion model to be usedfor a CGBF. The motion model can be developed as follows. Assume that atarget maneuvers at times t_(k), k=1, 2, 3, . . . These times areindependent of the measurement update times and are within an updateinterval. At each time t_(k), the target undergoes a constantacceleration and/or constant turn maneuver that persists until the nextmaneuver time. Assume, at time t_(k-1), the target is at location(x_(k-1), y_(k-1)) with speed and course, s_(k-1), θ_(k-1),respectively. Let {dot over (s)} be a (randomly) selected acceleration(along the direction of travel) and {dot over (θ)} be a (randomly)selected turn rate that the target will undergo for the duration of thenext movement. Let τ=t_(k)−t_(k-1) be the duration of the next maneuver.Then the target state as a function of time over the interval 0 to τwould be defined by the equations below:

$\begin{matrix}{{s(t)} = {s_{k - 1} + {\overset{.}{s}t}}} & {{Eq}.\mspace{14mu} 4} \\{{\theta (t)} = {\theta_{k - 1} + {\overset{.}{\theta}}_{t}}} & {{Eq}.\mspace{14mu} 5} \\{{\overset{.}{x}(t)} = {{s(t)}\sin \; {\theta (t)}}} & {{Eq}.\mspace{14mu} 6} \\{{\overset{.}{y}(t)} = {{s(t)}\cos \; {\theta (t)}}} & {{Eq}.\mspace{14mu} 7} \\\begin{matrix}{{x(t)} = {x_{k - 1} + {\int_{0}^{t}{{\overset{.}{x}(t)}\ {t}}}}} \\{= {{x_{k - 1} + {\frac{1}{{\overset{.}{\theta}}^{2}}\left\lbrack {{\overset{.}{s}\; \sin \; {\theta (t)}} - {{s(t)}\overset{.}{\theta}\; \cos \; {\theta (t)}}} \right\rbrack}}_{0}^{t}}} \\{= {x_{k - 1} + {\frac{1}{{\overset{.}{\theta}}^{2}}\begin{bmatrix}{{\overset{.}{s}\; \sin \; {\theta (t)}} - {{s(t)}\overset{.}{\theta}\; \cos \; \theta (t)}\; -} \\{{\overset{.}{s}\; \sin \; \theta_{k - 1}} + {s_{k - 1}\overset{.}{\theta}\; \cos \; \theta_{k - 1}}}\end{bmatrix}}}}\end{matrix} & {{Eq}.\mspace{14mu} 8} \\\begin{matrix}{{y(t)} = {y_{k - 1} + {\int_{0}^{t}{{\overset{.}{y}(t)}\ {t}}}}} \\{= {{y_{k - 1} + {\frac{1}{{\overset{.}{\theta}}^{2}}\left\lbrack {{\overset{.}{s}\; \cos \; {\theta (t)}} + {{s(t)}\overset{.}{\theta}\; \sin \; {\theta (t)}}} \right\rbrack}}_{0}^{t}}} \\{= {y_{k - 1} + {\frac{1}{{\overset{.}{\theta}}^{2}}\begin{bmatrix}{{\overset{.}{s}\; \cos \; {\theta (t)}} + {{s(t)}\overset{.}{\theta}\; \sin \; {\theta (t)}}\; -} \\{{\overset{.}{s}\; \cos \; \theta_{k - 1}} - {s_{k - 1}\overset{.}{\theta}\; \sin \; \theta_{k - 1}}}\end{bmatrix}}}}\end{matrix} & {{Eq}.\mspace{14mu} 9}\end{matrix}$

Even though the motion equations are mathematically sound, there couldbe computational problems when the maximum turn rate is (nearly) zero.To avoid these problems, let ε_(θ)>0 be the smallest turn rate value atwhich smaller values are equivalent to the larger moving along aconstant straight line course. Then the equations 8 and 9 become.

$\begin{matrix}{{x(t)} = {x_{k - 1} + \left\{ \begin{matrix}{\frac{1}{{\overset{.}{\theta}}^{2}}\begin{bmatrix}{{\overset{.}{s}\; \sin \; {\theta (t)}} - {{s(t)}\overset{.}{\theta}\; \cos \; \theta (t)}\; -} \\{{\overset{.}{s}\; \sin \; \theta_{k - 1}} + {s_{k - 1}\overset{.}{\theta}\; \cos \; \theta_{k - 1}}}\end{bmatrix}} & {{{if}\mspace{14mu} {\overset{.}{\theta}}} > \varepsilon_{\overset{.}{\theta}}} \\{\left( {{s_{k - 1}t} + {\frac{1}{2}\overset{.}{s}t^{2}}} \right){\sin \left( {\theta_{k - 1} + {\overset{.}{\theta}}_{\tau}} \right)}} & {otherwise}\end{matrix} \right.}} & {{Eq}.\mspace{14mu} 10} \\{{y(t)} = {y_{k - 1} + \left\{ \begin{matrix}{\frac{1}{{\overset{.}{\theta}}^{2}}\begin{bmatrix}{{\overset{.}{s}\; \cos \; {\theta (t)}} + {{s(t)}\overset{.}{\theta}\; \sin \; {\theta (t)}}\; -} \\{{\overset{.}{s}\; \cos \; \theta_{k - 1}} - {s_{k - 1}\overset{.}{\theta}\; \sin \; \theta_{k - 1}}}\end{bmatrix}} & {{{if}\mspace{14mu} {\overset{.}{\theta}}} > \varepsilon_{\overset{.}{\theta}}} \\{\left( {{s_{k - 1}t} + {\frac{1}{2}\overset{.}{s}t^{2}}} \right){\cos \left( {\theta_{k - 1} + {\overset{.}{\theta}}_{\tau}} \right)}} & {otherwise}\end{matrix} \right.}} & {{Eq}.\mspace{14mu} 11}\end{matrix}$

The Mean and Variance of the Transition Density (Single Maneuver)

Assume there is only one maneuver over the update time interval, i.e.,T=τ. Thus, the maneuver duration, τ, will be used to emphasize thisfact. Refer to equations 4-9. The state at the end of the maneuver isfound using t=τ. Recall that the subscript k refers to the (true) stateof the target after the k^(th) maneuver within an update time interval.A subscript of zero is the true state of the target at the beginning ofthe time interval. So:

$\begin{matrix}{\mspace{79mu} {s_{1} = {{s(\tau)} = {s_{0} + {\overset{.}{s}\tau}}}}} & {{Eq}.\mspace{14mu} 12} \\{\mspace{76mu} {\theta_{1} = {{\theta (\tau)} = {\theta_{0} + {\overset{.}{\theta}\tau}}}}} & {{Eq}.\mspace{14mu} 13} \\{\mspace{76mu} {{\overset{.}{x}}_{1} = {{\overset{.}{x}(\tau)} = {{{s(\tau)}\sin \; {\theta (\tau)}} = {s_{1}\sin \; \theta_{1}}}}}} & {{Eq}.\mspace{14mu} 14} \\{\mspace{76mu} {{\overset{.}{y}}_{1} = {{\overset{.}{y}(\tau)} = {{{s(\tau)}\cos \; {\theta (\tau)}} = {s_{1}\cos \; \theta_{1}}}}}} & {{Eq}.\mspace{14mu} 15} \\{x_{1} = {{x(\tau)} = {x_{0} + {\frac{1}{{\overset{.}{\theta}}^{2}}\left\lbrack {{\overset{.}{s}\sin \; \theta_{1}} - {s_{1}\overset{.}{\theta}\cos \; \theta_{1}} - {\overset{.}{s}\; \sin \; \theta_{0}} + {s_{0}\overset{.}{\theta}\cos \; \theta_{0}}} \right\rbrack}}}} & {{Eq}.\mspace{14mu} 16} \\{y_{1} = {{y(\tau)} = {y_{0} + {\frac{1}{{\overset{.}{\theta}}^{2}}\left\lbrack {{\overset{.}{s}\cos \; \theta_{1}} + {s_{1}\overset{.}{\theta}\sin \; \theta_{1}} - {\overset{.}{s}\; \cos \; \theta_{0}} - {s_{0}\overset{.}{\theta}\sin \; \theta_{0}}} \right\rbrack}}}} & {{Eq}.\mspace{14mu} 17}\end{matrix}$

The Mean Variance of the Speed and Course

The mean and variance of the speed and course is determined using theequations below:

s ₁ =s(τ)=s ₀ +{dot over (s)}τ and θ₁=θ(τ)=θ₀+{dot over (θ)}τ  Eq. 18

E[s ₁ ]=E[s ₀ +{dot over (s)}τ]=E[s ₀ ]+E[{dot over (s)}τ]=E[s ₀]+τE[{dot over (s)}]=s ₀   Eq. 19

And similarly from equations 13,

E[θ ₁ ]=E[θ _(o)+{dot over (θ)}τ]=E[θ ₀ ]+τE[{dot over (θ)}]=θ ₀   Eq.20

The expected speed and course after the target maneuvers are identicalto what they were prior to the maneuver. Calculation of the variancesfor the speed and course are calculated as follows:

Var[s ₁]=Var[s ₀+{dot over (s)}τ]=Var[s₀]+Var[{dot over (s)}τ]=Var[s₀]+τ²Var[{dot over (s)}]=τ²σ_(s) ²   Eq. 21

Var[θ₁]=Var[θ₀+{dot over (θ)}τ]=Var[θ₀]+τ²Var[{dot over (θ)}]=τ²σ_(θ) ²  Eq. 22

The selection of accelerations and turn rates were picked from a uniformdistribution over the range of possible rates. As a result,

${\sigma_{\overset{.}{s}}^{2} = {\frac{\left( {2\; A} \right)^{2}}{12} = {{\frac{A^{2}}{3}\mspace{14mu} {and}\mspace{14mu} \sigma_{\overset{.}{\theta}}^{2}} = {\frac{\left( {2\overset{.}{\Theta}} \right)^{2}}{12} = \frac{{\overset{.}{\Theta}}^{2}}{3}}}}},$

where A is the maximum acceleration and {dot over (Θ)} is the maximumturn rate. The variance of the velocity is quadratic with respect totime, or equivalently, the standard deviation of the velocity is linearin time.

Mean and Variance of the Component Velocities

The mean and variance of the component velocities can also bedetermined. In doing so, assume the distribution is in polar space andthe acceleration change is independent of the course change. Therefore,the following calculations are used.

E[{dot over (x)} ₁ ]=E[s ₁ sin θ₁]  Eq. 23

E[{dot over (y)} ₁ ]=E[s ₁ cos θ₁]  Eq. 24

Since it is assumed that the acceleration change is independent of thecourse change,

E[{dot over (x)} ₁ ]=E[s ₁ sin θ₁ ]=E[s ₁ ]E[sin θ₁]  Eq. 25

E[{dot over (y)} ₁ ]=E[s ₁ cos θ₁ ]=E[s ₁ ]E[cos θ₁]  Eq. 26

the E[s₁] was already found in equation 19. Using the properties of sineand cosine,

$\begin{matrix}\begin{matrix}{{E\left\lbrack {\sin \; \theta_{1}} \right\rbrack} = {E\left\lbrack {\sin\left( {\theta_{0} + {\overset{.}{\theta}\tau}} \right)} \right\rbrack}} \\{= {E\left\lbrack {{\sin \; \theta_{0}\cos \; \overset{.}{\theta}\tau} + {\sin \; \overset{.}{\theta}{\tau cos}\; \theta_{0}}} \right\rbrack}} \\{= {{\sin \; \theta_{0}{E\left\lbrack {\cos \; \overset{.}{\theta}\tau} \right\rbrack}} + {\cos \; \theta_{0}{E\left\lbrack {\sin \; \overset{.}{\theta}\tau} \right\rbrack}}}}\end{matrix} & {{Eq}.\mspace{14mu} 27} \\\begin{matrix}{{E\left\lbrack {\cos \; \theta_{1}} \right\rbrack} = {E\left\lbrack {\cos\left( {\theta_{0} + {\overset{.}{\theta}\tau}} \right)} \right\rbrack}} \\{= {E\left\lbrack {{\cos \; \theta_{0}\cos \; \overset{.}{\theta}\tau} - {\sin \; \theta_{0}\sin \; \overset{.}{\theta}\tau}} \right\rbrack}} \\{= {{\cos \; \theta_{0}{E\left\lbrack {\cos \; \overset{.}{\theta}\tau} \right\rbrack}} - {\sin \; \theta_{0}{E\left\lbrack {\sin \; \overset{.}{\theta}\tau} \right\rbrack}}}}\end{matrix} & {{Eq}.\mspace{14mu} 28}\end{matrix}$

The remaining expected values can be found using the definition ofexpected value. Assume, as before, that {dot over (θ)} is a uniformrandom variable on the interval [−{dot over (Θ)}, +{dot over (Θ)}].Therefore,

$\begin{matrix}\begin{matrix}{{E\left\lbrack {\sin \; \overset{.}{\theta}\tau} \right\rbrack} = {\int_{- \overset{.}{\Theta}}^{+ \overset{.}{\Theta}}{\left( \frac{1}{2\; \overset{.}{\Theta}} \right)\sin \; \overset{.}{\theta}\tau {\overset{.}{\theta}}}}} \\{= {{{- \left( \frac{1}{2\overset{.}{\Theta}\tau} \right)}\cos \; \overset{.}{\theta}\tau}|_{- \overset{.}{\Theta}}^{+ \overset{.}{\Theta}}}} \\{= {{- \left( \frac{1}{2\overset{.}{\Theta\tau}} \right)}\left( {{\cos \; \overset{.}{\Theta}\tau} - {\cos \; \overset{.}{\Theta}\tau}} \right)}} \\{= 0}\end{matrix} & {{Eq}.\mspace{14mu} 29} \\{\begin{matrix}{{E\left\lbrack {\cos \; \overset{.}{\theta}\tau} \right\rbrack} = {\int_{- \overset{.}{\Theta}}^{+ \overset{.}{\Theta}}{\left( \frac{1}{2\; \overset{.}{\Theta}} \right)\cos \; \overset{.}{\theta}\tau {\overset{.}{\theta}}}}} \\{= {{\left( \frac{1}{2\overset{.}{\Theta}\tau} \right)\sin \; \overset{.}{\theta}\tau}|_{- \overset{.}{\Theta}}^{+ \overset{.}{\Theta}}}} \\{= {\left( \frac{1}{2\overset{.}{\Theta\tau}} \right)\left( {{\sin \; \overset{.}{\Theta}\tau} + {\sin \; \overset{.}{\Theta}\tau}} \right)}} \\{= \frac{\sin \; \overset{.}{\Theta}\tau}{\overset{.}{\Theta}\tau}}\end{matrix}{{Thus},}} & {{Eq}.\mspace{14mu} 30} \\{{E\left\lbrack {\sin \; \theta_{1}} \right\rbrack} = {\sin \; \theta_{0}\frac{\sin \; \overset{.}{\Theta}\tau}{\overset{.}{\Theta}\tau}}} & {{Eq}.\mspace{14mu} 31} \\{{E\left\lbrack {\cos \; \theta_{1}} \right\rbrack} = {\cos \; \theta_{0}\frac{\sin \; \overset{.}{\Theta}\tau}{\overset{.}{\Theta}\tau}}} & {{Eq}.\mspace{14mu} 32}\end{matrix}$

From equations 27, 28, 29, 31, and 32,

$\begin{matrix}{{E\left\lbrack {\overset{.}{x}}_{1} \right\rbrack} = {s_{0}\sin \; \theta_{0}\frac{\sin \; \overset{.}{\Theta}\tau}{\overset{.}{\Theta}\tau}}} & {{Eq}.\mspace{14mu} 33} \\{{E\left\lbrack {\overset{.}{y}}_{1} \right\rbrack} = {s_{0}\cos \; \theta_{0}\frac{\sin \; \overset{.}{\Theta}\tau}{\overset{.}{\Theta}\tau}}} & {{Eq}.\mspace{14mu} 34}\end{matrix}$

Where {dot over (Θ)}τ is small,

$\frac{\sin \; \overset{.}{\Theta}\tau}{\overset{.}{\Theta}\tau} \approx 1.$

The variances of the component velocities can now be found.

Var[X]=E[X²]−(E[X])²   Eq. 35

Thus, from equations 27, 28, 33, 34,

$\begin{matrix}{{{Var}\left\lbrack {\overset{.}{x}}_{1} \right\rbrack} = {{{E\left\lbrack {\overset{.}{x}}_{1}^{2} \right\rbrack} - \left( {E\left\lbrack {\overset{.}{x}}_{1} \right\rbrack} \right)^{2}} = {{{E\left\lbrack s_{1}^{2} \right\rbrack}{E\left\lbrack {\sin^{2}\theta_{1}} \right\rbrack}} - \left( {s_{0}\sin \; \theta_{0}\frac{\sin \; \overset{.}{\Theta}\tau}{\overset{.}{\Theta}\tau}} \right)^{2}}}} & {{Eq}.\mspace{14mu} 36} \\{{{Var}\left\lbrack {\overset{.}{y}}_{1} \right\rbrack} = {{{E\left\lbrack {\overset{.}{y}}_{1}^{2} \right\rbrack} - \left( {E\left\lbrack {\overset{.}{y}}_{1} \right\rbrack} \right)^{2}} = {{{E\left\lbrack s_{1}^{2} \right\rbrack}{E\left\lbrack {\cos^{2}\theta_{1}} \right\rbrack}} - \left( {s_{0}\cos \; \theta_{0}\frac{\sin \; \overset{.}{\Theta}\tau}{\overset{.}{\Theta}\tau}} \right)^{2}}}} & {{Eq}.\mspace{14mu} 37}\end{matrix}$

Using equation E[s₁ ²]=Var[s₁]+(E[s₁])², so from equation 19 and 21,

E[s ₁ ²]=τ²σ_(s) ² +s ₀ ²   Eq. 38

E[sin² θ₁], and E[cos² θ₁] can be found using the definition of expectedvalue and Wolfram Mathematica Online Integrator.

$\begin{matrix}\begin{matrix}{{E\left\lbrack {\sin^{2}\theta_{1}} \right\rbrack} = {\int_{- \overset{.}{\Theta}}^{+ \overset{.}{\Theta}}{\left( \frac{1}{2\; \overset{.}{\Theta}} \right){\sin^{2}\left( {\theta_{0} + {\overset{.}{\theta}\tau}} \right)}{\overset{.}{\theta}}}}} \\{= {\frac{{2\left( {\theta_{0} + {\overset{.}{\theta}\tau}} \right)} - {\sin \; 2\left( {\theta_{0} + {\overset{.}{\theta}\tau}} \right)}}{8\; \overset{.}{\Theta}\tau}|_{- \overset{.}{\Theta}}^{+ \overset{.}{\Theta}}}} \\{= \frac{\begin{matrix}{{2\left( {\theta_{0} + {\overset{.}{\Theta}\tau}} \right)} - {\sin \; 2\left( {\theta_{0} + {\overset{.}{\Theta}\tau}} \right)} - {2\left( {\theta_{0} - {\overset{.}{\Theta}\tau}} \right)} +} \\{\sin \; 2\left( {\theta_{0} - {\overset{.}{\Theta}\tau}} \right)}\end{matrix}}{8\; \overset{.}{\Theta}\tau}} \\{= {\frac{1}{8\; \overset{.}{\Theta}\tau}\left\lbrack {{4\overset{.}{\Theta}\tau} - {\sin \; 2\left( {\theta_{0} + {\overset{.}{\Theta}\tau}} \right)} + {\sin \; 2\left( {\theta_{0} - {\overset{.}{\Theta}\tau}} \right)}} \right\rbrack}} \\{= {\frac{1}{2} + \frac{\begin{matrix}{{{- \sin}\; 2\theta_{0}\cos \; 2\; \overset{.}{\Theta}\tau} - {\cos \; 2\; \theta_{0}\sin \; 2\; \overset{.}{\Theta}\tau} +} \\{{\sin \; 2\; \theta_{0}\cos \; 2\overset{.}{\Theta}\tau} - {\cos \; 2\; \theta_{0}\sin \; 2\; \overset{.}{\Theta}\tau}}\end{matrix}}{8\; \overset{.}{\Theta}\tau}}} \\{= {\frac{1}{2} - \frac{\cos \; 2\; \theta_{0}\sin \; 2\overset{.}{\Theta}\tau}{4\; \overset{.}{\Theta\tau}}}} \\{= {\frac{1}{2}\left( {1 - {\frac{\sin \; 2\; \overset{.}{\Theta}\tau}{2\overset{.}{\Theta}\tau}\cos \; 2\; \theta_{0}}} \right)}}\end{matrix} & {{Eq}.\mspace{14mu} 39} \\\begin{matrix}{{E\left\lbrack {\cos^{2}\theta_{1}} \right\rbrack} = {\int_{- \overset{.}{\Theta}}^{+ \overset{.}{\Theta}}{\left( \frac{1}{2\; \overset{.}{\Theta}} \right){\cos^{2}\left( {\theta_{0} + {\overset{.}{\theta}\tau}} \right)}{\overset{.}{\theta}}}}} \\{= {\frac{{2\left( {\theta_{0} + {\overset{.}{\theta}\tau}} \right)} + {\sin \; 2\left( {\theta_{0} + {\overset{.}{\theta}\tau}} \right)}}{8\; \overset{.}{\Theta}\tau}|_{- \overset{.}{\Theta}}^{+ \overset{.}{\Theta}}}} \\{= \frac{\begin{bmatrix}{{2\left( {\theta_{0} + {\overset{.}{\Theta}\tau}} \right)} + {\sin \; 2\left( {\theta_{0} + {\overset{.}{\Theta}\tau}} \right)} - {2\left( {\theta_{0} - {\overset{.}{\Theta}\tau}} \right)} -} \\{\sin \; 2\left( {\theta_{0} - {\overset{.}{\Theta}\tau}} \right)}\end{bmatrix}}{8\; \overset{.}{\Theta}\tau}} \\{= {\frac{1}{8\; \overset{.}{\Theta}\tau}\left\lbrack {{4\overset{.}{\Theta}\tau} + {\sin \; 2\left( {\theta_{0} + {\overset{.}{\Theta}\tau}} \right)} - {\sin \; 2\left( {\theta_{0} - {\overset{.}{\Theta}\tau}} \right)}} \right\rbrack}} \\{= {\frac{1}{2} + \frac{\begin{matrix}{{\sin \; 2\theta_{0}\cos \; 2\; \overset{.}{\Theta}\tau} + {\cos \; 2\; \theta_{0}\sin \; 2\; \overset{.}{\Theta}\tau} -} \\{{\sin \; 2\; \theta_{0}\cos \; 2\overset{.}{\Theta}\tau} + {\cos \; 2\; \theta_{0}\sin \; 2\; \overset{.}{\Theta}\tau}}\end{matrix}}{8\; \overset{.}{\Theta}\tau}}} \\{= {\frac{1}{2} + \frac{\cos \; 2\; \theta_{0}\sin \; 2\overset{.}{\Theta}\tau}{4\; \overset{.}{\Theta\tau}}}} \\{= {\frac{1}{2}\left( {1 + {\frac{\sin \; 2\; \overset{.}{\Theta}\tau}{2\overset{.}{\Theta}\tau}\cos \; 2\; \theta_{0}}} \right)}}\end{matrix} & {{Eq}.\mspace{14mu} 40}\end{matrix}$

Therefore, using equations 36-40,

$\begin{matrix}{{{Var}\left\lbrack {\overset{.}{x}}_{1} \right\rbrack} = {{\left( {{\tau^{2}\sigma_{\overset{.}{s}}^{2}} + s_{0}^{2}} \right)\left( {\frac{1}{2} - \frac{\cos \; 2\; \theta_{0}\sin \; 2\; \overset{.}{\Theta}\tau}{4\overset{.}{\Theta}\tau}} \right)} - \left( {s_{0}\sin \; \theta_{0}\frac{\sin \; \overset{.}{\Theta}\tau}{\overset{.}{\Theta}\tau}} \right)^{2}}} & {{Eq}.\mspace{14mu} 41} \\{{{Var}\left\lbrack {\overset{.}{y}}_{1} \right\rbrack} = {{\left( {{\tau^{2}\sigma_{\overset{.}{s}}^{2}} + s_{0}^{2}} \right)\left( {\frac{1}{2} + \frac{\cos \; 2\; \theta_{0}\sin \; 2\; \overset{.}{\Theta}\tau}{4\overset{.}{\Theta}\tau}} \right)} - \left( {s_{0}\cos \; \theta_{0}\frac{\sin \; \overset{.}{\Theta}\tau}{\overset{.}{\Theta}\tau}} \right)^{2}}} & {{Eq}.\mspace{14mu} 42}\end{matrix}$

These equations emphasize how different the CGBF prediction model isfrom the one in the Kalman filter. To gain some appreciation for theseequations, consider the case when θ₀=0, i.e., the target's course is duenorth. Then, (41) and (42) become:

$\begin{matrix}{{{Var}\left\lbrack {\overset{.}{x}}_{1} \right\rbrack} = {\left( {{\tau^{2}\sigma_{\overset{.}{s}}^{2}} + s_{0}^{2}} \right)\left( {\frac{1}{2} - \frac{\sin \; 2\overset{.}{\Theta}\tau}{4\; \overset{.}{\Theta\tau}}} \right)}} & {{Eq}.\mspace{14mu} 43} \\{{{Var}\left\lbrack {\overset{.}{y}}_{1} \right\rbrack} = {{\left( {{\tau^{2}\sigma_{\overset{.}{s}}^{2}} + s_{0}^{2}} \right)\left( {\frac{1}{2} + \frac{\sin \; 2\overset{.}{\Theta}\tau}{4\; \overset{.}{\Theta\tau}}} \right)} - \left( {s_{0}\frac{\sin \; \overset{.}{\Theta}\tau}{\overset{.}{\Theta}\tau}} \right)^{2}}} & {{Eq}.\mspace{14mu} 44}\end{matrix}$

As {dot over (Θ)}τ goes to zero, i.e., {dot over (Θ)}τ→0, the ration

${\frac{\sin \; \overset{.}{\Theta}\tau}{\overset{.}{\Theta}\tau}->1},$

so:

$\begin{matrix}{{{{Var}\left\lbrack {\overset{.}{x}}_{1} \right\rbrack}->{\left( {{\tau^{2}\sigma_{\overset{.}{s}}^{2}} + s_{0}^{2}} \right)\left( {\frac{1}{2} - \frac{1}{2}} \right)}} = 0} & {{Eq}.\mspace{14mu} 45} \\{{{{Var}\left\lbrack {\overset{.}{y}}_{1} \right\rbrack}->{{\left( {{\tau^{2}\sigma_{\overset{.}{s}}^{2}} + s_{0}^{2}} \right)\left( {\frac{1}{2} + \frac{1}{2}} \right)} - \left( s_{0} \right)^{2}}} = {\tau^{2}\sigma_{\overset{.}{s}}^{2}}} & {{Eq}.\mspace{14mu} 46}\end{matrix}$

The variance in the x direction should go to zero because when {dot over(Θ)}τ=0, the target is (only) moving exactly along its initial course,so there is no uncertainty in course. The resulting variance in the ydirection (i.e., the variance in position uncertainty along the target'sinitial course) agrees with the process model for the Kalman filter.Thus, the CGBF uses a directional process model.

The Mean and Variance of the Position

The mean and variance of the position can also be determined. The meanof the position is considered first. Without any loss of generality,assume the target is initially moving due north. This assumption doesnot limit the generality because the distribution is being computedalong the initial course of the target. Thus, any other initial courseis just a rotation of the mean position to align with the target'sinitial course. With this assumption,

θ₀=0.

The distribution must be symmetric about the initial course since thetarget can turn left just as easily as turning right. Since the targetis assumed to have an initial course of zero, i.e., due north, thedistribution must be symmetric about the y-axis. Therefore, E[x₁]=x₀.However, this result will be proven. The E[x₁] is computed from themodified position Eq. 10 and 11.

$\begin{matrix}{{E\left\lbrack x_{1} \right\rbrack} = {E\left\lbrack {x_{0} + \left\{ \begin{matrix}{\frac{1}{{\overset{.}{\theta}}^{2}}\left\lbrack {{\overset{.}{s}\sin \; \theta_{1}} - {s_{1}\overset{.}{\theta}\cos \; \theta_{1}} - {\overset{.}{s}\sin \; \theta_{0}} + {s_{0}\overset{.}{\theta}\cos \; \theta_{0}}} \right\rbrack} & {{{if}\mspace{14mu} {\overset{.}{\theta}}} > \varepsilon_{\overset{.}{\theta}}} \\{\left( {{s_{0}\tau} + {\frac{1}{2}\overset{.}{s}\tau^{2}}} \right){\sin \left( {\theta_{0} + {\overset{.}{\theta}\tau}} \right)}} & {otherwise}\end{matrix} \right\rbrack} \right.}} & {{Eq}.\mspace{14mu} 47} \\{{E\left\lbrack y_{1} \right\rbrack} = {E\left\lbrack {y_{0} + \left\{ \begin{matrix}{\frac{1}{{\overset{.}{\theta}}^{2}}\left\lbrack {{\overset{.}{s}\cos \; \theta_{1}} + {s_{1}\overset{.}{\theta}\sin \; \theta_{1}} - {\overset{.}{s}\cos \; \theta_{0}} - {s_{0}\overset{.}{\theta}\sin \; \theta_{0}}} \right\rbrack} & {{{if}\mspace{14mu} {\overset{.}{\theta}}} > \varepsilon_{\overset{.}{\theta}}} \\{\left( {{s_{0}\tau} + {\frac{1}{2}\overset{.}{s}\tau^{2}}} \right){\cos \left( {\theta_{0} + {\overset{.}{\theta}\tau}} \right)}} & {otherwise}\end{matrix} \right\rbrack} \right.}} & {{Eq}.\mspace{14mu} 48}\end{matrix}$

As before, assume the target is initially moving due north, so θ₀=0.With this simplification:

$\begin{matrix}{{E\left\lbrack x_{1} \right\rbrack} = {E\left\lbrack {x_{0} + \left\{ \begin{matrix}{\frac{1}{{\overset{.}{\theta}}^{2}}\left\lbrack {{\overset{.}{s}\sin \; \overset{.}{\theta}\tau} - {s_{1}\overset{.}{\theta}\cos \; \overset{.}{\theta}\tau} + {s_{0}\overset{.}{\theta}}} \right\rbrack} & {{{if}\mspace{14mu} {\overset{.}{\theta}}} > \varepsilon_{\overset{.}{\theta}}} \\{\left( {{s_{0}\tau} + {\frac{1}{2}\overset{.}{s}\tau^{2}}} \right)\sin \; \overset{.}{\theta}\tau} & {otherwise}\end{matrix} \right\rbrack} \right.}} & {{Eq}.\mspace{14mu} 49} \\{{E\left\lbrack y_{1} \right\rbrack} = {E\left\lbrack {y_{0} + \left\{ \begin{matrix}{\frac{1}{{\overset{.}{\theta}}^{2}}\left\lbrack {{\overset{.}{s}\cos \; \overset{.}{\theta}\tau} + {s_{1}\overset{.}{\theta}\sin \; \overset{.}{\theta}\tau} - \overset{.}{s}} \right\rbrack} & {{{if}\mspace{14mu} {\overset{.}{\theta}}} > \varepsilon_{\overset{.}{\theta}}} \\{\left( {{s_{0}\tau} + {\frac{1}{2}\overset{.}{s}\tau^{2}}} \right)\cos \; \overset{.}{\theta}\tau} & {otherwise}\end{matrix} \right\rbrack} \right.}} & {{Eq}.\mspace{14mu} 50}\end{matrix}$

And again assuming {dot over (s)} and {dot over (θ)} are independentrandom variables,

$\begin{matrix}{{E\left\lbrack x_{1} \right\rbrack} = {x_{0} + \left\{ \begin{matrix}{{{E\left\lbrack \overset{.}{s} \right\rbrack}{E\left\lbrack \frac{\sin \; \overset{.}{\theta}\tau}{{\overset{.}{\theta}}^{2}} \right\rbrack}} - {{E\left\lbrack s_{1} \right\rbrack}{E\left\lbrack \frac{\cos \; \overset{.}{\theta}\tau}{\overset{.}{\theta}} \right\rbrack}} + {s_{0}{E\left\lbrack \frac{1}{\overset{.}{\theta}} \right\rbrack}}} & {{{if}\mspace{14mu} {\overset{.}{\theta}}} > \varepsilon_{\overset{.}{\theta}}} \\{{E\left\lbrack {{s_{0}\tau} + {\frac{1}{2}\overset{.}{s}\tau^{2}}} \right\rbrack}{E\left\lbrack {\sin \; \overset{.}{\theta}\tau} \right\rbrack}} & {otherwise}\end{matrix} \right.}} & {{Eq}.\mspace{14mu} 51} \\{{E\left\lbrack y_{1} \right\rbrack} = {y_{0} + \left\{ \begin{matrix}{{{E\left\lbrack \overset{.}{s} \right\rbrack}{E\left\lbrack \frac{\cos \; \overset{.}{\theta}\tau}{{\overset{.}{\theta}}^{2}} \right\rbrack}} + {{E\left\lbrack s_{1} \right\rbrack}{E\left\lbrack \frac{\sin \; \overset{.}{\theta}\tau}{\overset{.}{\theta}} \right\rbrack}} - {{E\left\lbrack \overset{.}{s} \right\rbrack}{E\left\lbrack \frac{1}{{\overset{.}{\theta}}^{2}} \right\rbrack}}} & {{{if}\mspace{14mu} {\overset{.}{\theta}}} > \varepsilon_{\overset{.}{\theta}}} \\{{E\left\lbrack {{s_{0}\tau} + {\frac{1}{2}\overset{.}{s}\tau^{2}}} \right\rbrack}{E\left\lbrack {\cos \; \overset{.}{\theta}\tau} \right\rbrack}} & {otherwise}\end{matrix} \right.}} & {{Eq}.\mspace{14mu} 52}\end{matrix}$

Since E[{dot over (s)}]=0,

$\begin{matrix}{{E\left\lbrack x_{1} \right\rbrack} = {x_{0} + \left\{ {\begin{matrix}{{{- s_{0}}{E\left\lbrack \frac{\cos \; \overset{.}{\theta}\tau}{\overset{.}{\theta}} \right\rbrack}} + {s_{0}{E\left\lbrack \frac{1}{\overset{.}{\theta}} \right\rbrack}}} & {{{if}{\overset{.}{\theta}}} > \varepsilon_{\overset{.}{\theta}}} \\{s_{0}\tau \; {E\left\lbrack {\sin \; \overset{.}{\theta}\tau} \right\rbrack}} & {otherwise}\end{matrix} = {x_{0} + \left\{ \begin{matrix}{s_{0}{E\left\lbrack \frac{1 - {\cos \; \overset{.}{\theta}\tau}}{\overset{.}{\theta}} \right\rbrack}} & {{{if}\mspace{14mu} {\overset{.}{\theta}}} > \varepsilon_{\overset{.}{\theta}}} \\{s_{0}\tau \; {E\left\lbrack {\sin \; \overset{.}{\theta}\tau} \right\rbrack}} & {otherwise}\end{matrix} \right.}} \right.}} & {{Eq}.\mspace{14mu} 53} \\{\mspace{79mu} {{E\left\lbrack y_{1} \right\rbrack} = {y_{0} + \left\{ \begin{matrix}{s_{0}{E\left\lbrack \frac{\sin \; \overset{.}{\theta}\tau}{\overset{.}{\theta}} \right\rbrack}} & {{{if}{\overset{.}{\theta}}} > \varepsilon_{\overset{.}{\theta}}} \\{s_{0}\tau \; {E\left\lbrack {\cos \; \overset{.}{\theta}\tau} \right\rbrack}} & {otherwise}\end{matrix} \right.}}} & {{Eq}.\mspace{14mu} 54}\end{matrix}$

The expected values in Eq. 51 and Eq. 52 can be written as:

$\begin{matrix}{{E\left\lbrack x_{1} \right\rbrack} = {x_{0} + {s_{0}\begin{bmatrix}{{\int_{- \overset{.}{\Theta}}^{- \varepsilon_{\overset{.}{\theta}}}{\left( \frac{1}{2\overset{.}{\Theta}} \right)\left( \frac{1 - {\cos \; \overset{.}{\theta}\tau}}{\overset{.}{\theta}} \right){\overset{.}{\theta}}}} +} \\{{\tau {\int_{- \varepsilon_{\overset{.}{\theta}}}^{+ \varepsilon_{\overset{.}{\theta}}}{\left( \frac{1}{2\; \overset{.}{\Theta}} \right)\sin \; \overset{.}{\theta}\tau {\overset{.}{\theta}}}}} +} \\{\int_{+ \varepsilon_{\overset{.}{\theta}}}^{+ \overset{.}{\Theta}}{\left( \frac{1}{2\overset{.}{\Theta}} \right)\left( \frac{1 - {\cos \; \overset{.}{\theta}\tau}}{\overset{.}{\theta}} \right){\overset{.}{\theta}}}}\end{bmatrix}}}} & {{Eq}.\mspace{14mu} 55} \\{{E\left\lbrack y_{1} \right\rbrack} = {y_{0} + {s_{0}\begin{bmatrix}{{\int_{- \overset{.}{\Theta}}^{- \varepsilon_{\overset{.}{\theta}}}{\left( \frac{\sin \; \overset{.}{\theta}\tau}{\overset{.}{\theta}} \right)\left( \frac{1}{2\; \overset{.}{\Theta}} \right){\overset{.}{\theta}}}} +} \\{{\int_{- \varepsilon_{\overset{.}{\theta}}}^{+ \varepsilon_{\overset{.}{\theta}}}{\left( \frac{1}{{2\overset{.}{\Theta}}\;} \right)\cos \; \overset{.}{\theta}\tau {\overset{.}{\theta}}}} +} \\{\tau {\int_{- \varepsilon_{\overset{.}{\theta}}}^{+ \overset{.}{\Theta}}{\left( \frac{\sin \; \overset{.}{\theta}\tau}{\overset{.}{\theta}} \right)\left( \frac{1}{2\overset{.}{\Theta}} \right){\overset{.}{\theta}}}}}\end{bmatrix}}}} & {{Eq}.\mspace{14mu} 56}\end{matrix}$

These integrals are identical to the expected values for E[sin {dot over(θ)}τ], E]cos {dot over (θ)}τ],

${E\left\lbrack \frac{\sin \; \overset{.}{\theta}\tau}{\overset{.}{\theta}} \right\rbrack},{and}$$E\left\lbrack \frac{1 - {\cos \; \overset{.}{\theta}\tau}}{\overset{.}{\theta}} \right\rbrack$

that were found earlier, except now with different limits ofintegration. Within the limits now defined, each function iswell-behaved with no infinities in the range. First consider E[x₁]. Allthree integrals in Eq. 55 do not need to be evaluated. Since theintegrands are odd functions and symmetrical area is being integrated,the area must sum to zero. Therefore, all the integrals for E[x₁] arezero. Now consider E[y₁]. The second integral in Eq. 56 is identical toE[cos {dot over (θ)}τ] that was already evaluated in (Eq. 30) exceptwith different limits:

$\begin{matrix}{{\int_{- \varepsilon_{\overset{.}{\theta}}}^{+ \varepsilon_{\overset{.}{\theta}}}{\left( \frac{1}{2\overset{.}{\Theta}} \right)\cos \; \overset{.}{\theta}\tau {\overset{.}{\theta}}}} = {{{\left( \frac{1}{2\overset{.}{\Theta}\tau} \right)\sin \; \overset{.}{\theta}\tau}|_{- \varepsilon_{\overset{.}{\theta}}}^{+ \varepsilon_{\overset{.}{\theta}}}} = {{\left( \frac{1}{2\overset{.}{\Theta}\tau} \right)\left( {{\sin \; \varepsilon_{\overset{.}{\theta}}\tau} + {\sin \; \varepsilon_{\overset{.}{\theta}}\tau}} \right)} = \frac{\sin \; \varepsilon_{\overset{.}{\theta}}\tau}{\overset{.}{\Theta}\tau}}}} & {{Eq}.\mspace{14mu} 57}\end{matrix}$

The first and last integrals in Eq. 56 are similar to evaluating E[sin{dot over (θ)}τ/{dot over (θ)}].

$\begin{matrix}{{{\int_{- \overset{.}{\Theta}}^{- \varepsilon_{\overset{.}{\theta}}}{\left( \frac{\sin \; \overset{.}{\theta}\tau}{\overset{.}{\theta}} \right)\left( \frac{1}{2\; \overset{.}{\Theta}} \right){\overset{.}{\theta}}}} + {\int_{+ \varepsilon_{\overset{.}{\theta}}}^{+ \overset{.}{\Theta}}{\left( \frac{\sin \; \overset{.}{\theta}\tau}{\overset{.}{\theta}} \right)\left( \frac{1}{2\overset{.}{\Theta}} \right){\overset{.}{\theta}}}}} = {{\frac{{Si}\left( {\overset{.}{\theta}\tau} \right)}{2\overset{.}{\Theta}}{_{- \overset{.}{\Theta}}^{- \varepsilon_{\overset{.}{\theta}}}{+ \frac{{Si}\left( {\overset{.}{\theta}\tau} \right)}{2\overset{.}{\Theta}}}}_{+ \varepsilon_{\overset{.}{\theta}}}^{+ \overset{.}{\Theta}}} = {{\frac{1}{2\overset{.}{\Theta}}\begin{bmatrix}{{{Si}\left( {{- \varepsilon_{\overset{.}{\theta}}}\tau} \right)} - {{Si}\left( {{- \overset{.}{\Theta}}\tau} \right)} +} \\{{{Si}\left( {\overset{.}{\Theta}\tau} \right)} - {{Si}\left( {\varepsilon_{\overset{.}{\theta}}\tau} \right)}}\end{bmatrix}} = {{\frac{1}{2\overset{.}{\Theta}}\begin{bmatrix}{{- {{Si}\left( {\varepsilon_{\overset{.}{\theta}}\tau} \right)}} + {{Si}\left( {\overset{.}{\Theta}\tau} \right)} +} \\{{{Si}\left( {\overset{.}{\Theta}\tau} \right)} - {{Si}\left( {\varepsilon_{\overset{.}{\theta}}\tau} \right)}}\end{bmatrix}} = {\frac{1}{\overset{.}{\Theta}}\left\lbrack {{{Si}\left( {\overset{.}{\Theta}\tau} \right)} - {{Si}\left( {\varepsilon_{\overset{.}{\theta}}\tau} \right)}} \right\rbrack}}}}} & {{Eq}.\mspace{14mu} 58}\end{matrix}$

Using these results in Eq. 55 and Eq. 56 the expected position using themodified motion equation is:

$\begin{matrix}{{E\left\lbrack x_{1} \right\rbrack} = x_{0}} & {{Eq}.\mspace{14mu} 59} \\{{E\left\lbrack y_{1} \right\rbrack} = {y_{0} + {\frac{S_{k - 1}}{\overset{.}{\Theta}}\left\lbrack {{{Si}\left( {\overset{.}{\Theta}\tau} \right)} - {{Si}\left( {\varepsilon_{\overset{.}{\theta}}\tau} \right)} + {\sin \; \varepsilon_{\overset{.}{\theta}}\tau}} \right\rbrack}}} & {{Eq}.\mspace{14mu} 60}\end{matrix}$

Since ∈_({dot over (θ)}) is small, both Si(∈_({dot over (θ)})τ)≈0 andsin ∈_({dot over (θ)})τ≈0, so

${E\left\lbrack y_{1} \right\rbrack} \approx {y_{0} + {s_{0}{\frac{{Si}\left( {\overset{.}{\Theta}\tau} \right)}{\overset{.}{\Theta}}.}}}$

Now consider the variance of the position. To avoid the issue ofpossible division by zero, the modification, Eq. 10 and 11 or 49 and 50,as well will be used to find the variance.

$\begin{matrix}{{E\left\lbrack x_{1}^{2} \right\rbrack} = {E\left\lbrack \left( {x_{0} + \left\{ \begin{matrix}{\frac{1}{{\overset{.}{\theta}}^{2}}\begin{bmatrix}{{\overset{.}{s}\sin \; \theta_{1}} - {s_{1}\overset{.}{\theta}\cos \; \theta_{1}} -} \\{{\overset{.}{s}\sin \; \theta_{0}} + {s_{0}\overset{.}{\theta}\cos \; \theta_{0}}}\end{bmatrix}} & {{{if}\mspace{14mu} {\overset{.}{\theta}}} > \varepsilon_{\overset{.}{\theta}}} \\{\left( {{s_{0}\tau} + {\frac{1}{2}\overset{.}{s}\tau^{2}}} \right){\sin \left( {\theta_{0} + {\overset{.}{\theta}\tau}} \right)}} & {otherwise}\end{matrix} \right)^{2}} \right\rbrack \right.}} & {{Eq}.\mspace{14mu} 61} \\{{E\left\lbrack y_{1}^{2} \right\rbrack} = {E\left\lbrack \left( {y_{0} + \left\{ \begin{matrix}{\frac{1}{{\overset{.}{\theta}}^{2}}\begin{bmatrix}{{\overset{.}{s}\cos \; \theta_{1}} - {s_{1}\overset{.}{\theta}\sin \; \theta_{1}} -} \\{{\overset{.}{s}\cos \; \theta_{0}} + {s_{0}\overset{.}{\theta}\sin \; \theta_{0}}}\end{bmatrix}} & {{{if}\mspace{14mu} {\overset{.}{\theta}}} > \varepsilon_{\overset{.}{\theta}}} \\{\left( {{s_{0}\tau} + {\frac{1}{2}\overset{.}{s}\tau^{2}}} \right){\sin \left( {\theta_{0} + {\overset{.}{\theta}\tau}} \right)}} & {otherwise}\end{matrix} \right)^{2}} \right\rbrack \right.}} & {{Eq}.\mspace{14mu} 62}\end{matrix}$

All of the integrals arising from these expected values are solvable butproduce “messy” solutions. Thus, although a closed form solution can befound for the variance of the position, it is much too complicated to beuseful. Therefore, the predicted state and covariance will be estimatedusing numerical methods as described next.

The Mean and Variance of the Transition Density (Double Maneuver)

Now that the distribution parameters have been determined for the casewhen there is one maneuver over the update time interval, the case whenthere are two maneuvers is considered next. For this two-step case,assume the target starts with some randomly selected acceleration andturn rate and executes that maneuver for τ time, just as was assumed forthe one-step case. But at time τ, the target randomly selects anotheracceleration and turn rate and executes that maneuver for the remainderof the update time interval, T−τ. This assumed process now makes τ arandom variable as well. Maintaining consistency with the rest of theanalysis, assume τ is a uniform random variable such that 0≦τ≦T.Therefore,

$\begin{matrix}{{E\lbrack\tau\rbrack} = {\frac{1}{2}T}} & {{Eq}.\mspace{14mu} 63} \\{{{Var}\lbrack\tau\rbrack} = {\sigma_{\tau}^{2} = {\frac{1}{12}T^{2}}}} & {{Eq}.\mspace{14mu} 64}\end{matrix}$

To deal with the multiple maneuvers, the (randomly) selectedaccelerations and turn rates need to be subscripted as well. As before,assume at a specified update time, the target is at (x₀, y₀). Let thisupdate time be at time zero. The (predicted) state of the target isdesired for the next update time, T. Let τ be the time when the targetfinishes its first maneuver and starts its second (and last) maneuverover the time update interval. As before, the subscripts reflect targetstates within an update time interval. So (x₀, y₀) is the position ofthe target at the begin inning of the update time interval havinginitial speed s₀ and initial course θ₀. An acceleration {dot over (s)}₁and turn rate, {dot over (θ)}₁, is (randomly) selected. The target wouldthen have position (x₁, y₁) τ time later with speed s₁ and course θ₁. Atthis point in time, the target (randomly) selects a new acceleration andturn rate, {dot over (s)}₂ and {dot over (θ)}₂, respectively. The targetexecutes that maneuver for the remainder of the update time interval,T−τ. At the end of the update time period, T, the target would haveposition (x₂, y₂) with speed and course s₂ and θ₂, respectively. Usingthis notation,

$\begin{matrix}{s_{1} = {s_{0} + {{\overset{.}{s}}_{1}\tau}}} & {{Eq}.\mspace{14mu} 65} \\{s_{2} = {{s_{1} + {{\overset{.}{s}}_{2}\left( {T - \tau} \right)}} = {s_{0} + {{\overset{.}{s}}_{1}\tau} + {{\overset{.}{s}}_{2}\left( {T - \tau} \right)}}}} & {{Eq}.\mspace{14mu} 66} \\{\theta_{1} = {\theta_{0} + {{\overset{.}{\theta}}_{1}\tau}}} & {{Eq}.\mspace{14mu} 67} \\{\theta_{2} = {{\theta_{1} + {{\overset{.}{\theta}}_{2}\left( {T - \tau} \right)}} = {\theta_{0} + {{\overset{.}{\theta}}_{1}\tau} + {{\overset{.}{\theta}}_{2}\left( {T - \tau} \right)}}}} & {{Eq}.\mspace{14mu} 68} \\{{\overset{.}{x}}_{1} = {s_{1}\sin \; \theta_{1}}} & {{Eq}.\mspace{14mu} 69} \\{{\overset{.}{x}}_{2} = {s_{2}\sin \; \theta_{2}}} & {{Eq}.\mspace{14mu} 70} \\{{\overset{.}{y}}_{1} = {s_{1}\cos \; \theta_{1}}} & {{Eq}.\mspace{14mu} 71} \\{{\overset{.}{y}}_{2} = {s_{2}\cos \; \theta_{2}}} & {{Eq}.\mspace{14mu} 72} \\{x_{1} = {x_{0} + {\frac{1}{{\overset{.}{\theta}}_{1}^{2}}\left\lbrack {{{\overset{.}{s}}_{1}\sin \; \theta_{1}} - {s_{1}{\overset{.}{\theta}}_{1}\cos \; \theta_{1}} - {{\overset{.}{s}}_{1}\sin \; \theta_{0}} + {s_{0}{\overset{.}{\theta}}_{1}\cos \; \theta_{0}}} \right\rbrack}}} & {{Eq}.\mspace{14mu} 73} \\{x_{2} = {x_{1} + {\frac{1}{{\overset{.}{\theta}}_{2}^{2}}\left\lbrack {{{\overset{.}{s}}_{2}\sin \; \theta_{2}} - {s_{2}{\overset{.}{\theta}}_{2}\cos \; \theta_{2}} - {{\overset{.}{s}}_{2}\sin \; \theta_{1}} + {s_{1}{\overset{.}{\theta}}_{2}\cos \; \theta_{1}}} \right\rbrack}}} & {{Eq}.\mspace{14mu} 74} \\{y_{1} = {y_{0} + {\frac{1}{{\overset{.}{\theta}}_{1}^{2}}\left\lbrack {{{\overset{.}{s}}_{1}\cos \; \theta_{1}} - {s_{1}{\overset{.}{\theta}}_{1}\sin \; \theta_{1}} - {{\overset{.}{s}}_{1}\cos \; \theta_{0}} + {s_{0}{\overset{.}{\theta}}_{1}\sin \; \theta_{0}}} \right\rbrack}}} & {{Eq}.\mspace{14mu} 75} \\{y_{2} = {y_{1} + {\frac{1}{{\overset{.}{\theta}}_{2}^{2}}\left\lbrack {{{\overset{.}{s}}_{2}\cos \; \theta_{2}} - {s_{2}{\overset{.}{\theta}}_{2}\sin \; \theta_{2}} - {{\overset{.}{s}}_{2}\cos \; \theta_{1}} + {s_{1}{\overset{.}{\theta}}_{2}\sin \; \theta_{1}}} \right\rbrack}}} & {{Eq}.\mspace{14mu} 76}\end{matrix}$

The Mean and Variance of the Speed and Course

To find the mean and variance of the speed and course for thetwo-maneuver case, note that it is necessary to deal with a product ofindependent random variables. It is known that the variance of a productof independent random variables, A, B is given by:

Var[AB]=(E[A])²Var[B]+(E[B])²Var[A]+Var[A]Var[B]  Eq. 77

Using this property along with

${{E\left\lbrack {\overset{.}{s}}_{1} \right\rbrack} = {{E\left\lbrack {\overset{.}{s}}_{2} \right\rbrack} = {{E\left\lbrack {\overset{.}{\theta}}_{1} \right\rbrack} = {{E\left\lbrack {\overset{.}{\theta}}_{2} \right\rbrack} = 0}}}},{{E\left\lbrack {T - \tau} \right\rbrack} = {{E\lbrack\tau\rbrack} = {\frac{1}{2}T}}},{{{Var}\left\lbrack {\overset{.}{s}}_{1} \right\rbrack} = {{{Var}\left\lbrack {\overset{.}{s}}_{2} \right\rbrack} = \sigma_{\overset{.}{s}}^{2}}},{and}$${{{Var}\left\lbrack {T - \tau} \right\rbrack} = {{{Var}\lbrack\tau\rbrack} = {\sigma_{\tau}^{2} = {\frac{1}{12}T^{2}}}}},{then}$$\begin{matrix}\begin{matrix}{{E\left\lbrack s_{2} \right\rbrack} = {{E\left\lbrack s_{0} \right\rbrack} + {E\left\lbrack {{\overset{.}{s}}_{1}\tau} \right\rbrack} + {E\left\lbrack {{\overset{.}{s}}_{2}\left( {T - \tau} \right)} \right\rbrack}}} \\{= {s_{0} + {{E\left\lbrack {\overset{.}{s}}_{1} \right\rbrack}{E\lbrack\tau\rbrack}} + {{E\left\lbrack {\overset{.}{s}}_{2} \right\rbrack}{E\left\lbrack {T - \tau} \right\rbrack}}}} \\{= s_{0}}\end{matrix} & {{Eq}.\mspace{14mu} 78} \\\begin{matrix}{{{Var}\left\lbrack s_{2} \right\rbrack} = {{{Var}\left\lbrack s_{0} \right\rbrack} + {{Var}\left\lbrack {{\overset{.}{s}}_{1}\tau} \right\rbrack} + {{Var}\left\lbrack {{\overset{.}{s}}_{2}\left( {T - \tau} \right)} \right\rbrack}}} \\{= {\begin{Bmatrix}{{\left( {E\left\lbrack {\overset{.}{s}}_{2} \right\rbrack} \right)^{2}{{Var}\lbrack\tau\rbrack}} + {\left( {E\lbrack\tau\rbrack} \right)^{2}{{Var}\left\lbrack {\overset{.}{s}}_{1} \right\rbrack}} +} \\{{{Var}\left\lbrack {\overset{.}{s}}_{1} \right\rbrack}{{Var}\lbrack\tau\rbrack}}\end{Bmatrix} +}} \\{\begin{Bmatrix}{{\left( {E\left\lbrack {\overset{.}{s}}_{2} \right\rbrack} \right)^{2}{{Var}\left\lbrack {T - \tau} \right\rbrack}} +} \\{{\left( {E\left\lbrack {T - \tau} \right\rbrack} \right)^{2}{{Var}\left\lbrack {\overset{.}{s}}_{2} \right\rbrack}} + {{{Var}\left\lbrack {\overset{.}{s}}_{2} \right\rbrack}{{Var}\left\lbrack {T - \tau} \right\rbrack}}}\end{Bmatrix}} \\{= {{2\left\lbrack {{\left( {\frac{1}{2}T} \right)^{2}\sigma_{\overset{.}{s}}^{2}} + {\sigma_{\overset{.}{s}}^{2}\left( {\frac{1}{12}T^{2}} \right)}}\; \right\rbrack} = {\frac{2}{3}T^{2}\sigma_{\overset{.}{s}}^{2}}}}\end{matrix} & {{Eq}.\mspace{14mu} 79} \\\begin{matrix}{{E\left\lbrack \theta_{2} \right\rbrack} = {{E\left\lbrack \theta_{0} \right\rbrack} + {E\left\lbrack {{\overset{.}{\theta}}_{1}\tau} \right\rbrack} + {E\left\lbrack {{\overset{.}{\theta}}_{2}\left( {T - \tau} \right)} \right\rbrack}}} \\{= {\theta_{0} + {{E\left\lbrack {\overset{.}{\theta}}_{1} \right\rbrack}{E\lbrack\tau\rbrack}} + {{E\left\lbrack {\overset{.}{\theta}}_{2} \right\rbrack}{E\left\lbrack {T - \tau} \right\rbrack}}}} \\{= \theta_{0}}\end{matrix} & {{Eq}.\mspace{14mu} 80} \\\begin{matrix}{{{Var}\left\lbrack \theta_{2} \right\rbrack} = {{{Var}\left\lbrack \theta_{0} \right\rbrack} + {{Var}\left\lbrack {{\overset{.}{\theta}}_{1}\tau} \right\rbrack} + {{Var}\left\lbrack {{\overset{.}{\theta}}_{2}\left( {T - \tau} \right)} \right\rbrack}}} \\{= {\begin{Bmatrix}{{\left( {E\left\lbrack {\overset{.}{\theta}}_{1} \right\rbrack} \right)^{2}{{Var}\lbrack\tau\rbrack}} + {\left( {E\lbrack\tau\rbrack} \right)^{2}{{Var}\left\lbrack {\overset{.}{\theta}}_{1} \right\rbrack}} +} \\{{{Var}\left\lbrack {\overset{.}{\theta}}_{1} \right\rbrack}{{Var}\lbrack\tau\rbrack}}\end{Bmatrix} +}} \\{\begin{Bmatrix}{{\left( {E\left\lbrack {\overset{.}{\theta}}_{2} \right\rbrack} \right)^{2}{{Var}\left\lbrack {T - \tau} \right\rbrack}} +} \\{{\left( {E\left\lbrack {T - \tau} \right\rbrack} \right)^{2}{{Var}\left\lbrack {\overset{.}{\theta}}_{2} \right\rbrack}} + {{{Var}\left\lbrack {\overset{.}{\theta}}_{2} \right\rbrack}{{Var}\left\lbrack {T - \tau} \right\rbrack}}}\end{Bmatrix}} \\{= {{2\left\lbrack {{\left( {\frac{1}{2}T} \right)^{2}\sigma_{\overset{.}{\theta}}^{2}} + {\sigma_{\overset{.}{\theta}}^{2}\left( {\frac{1}{12}T^{2}} \right)}}\; \right\rbrack} = {\frac{2}{3}T^{2}\sigma_{\overset{.}{\theta}}^{2}}}}\end{matrix} & {{Eq}.\mspace{14mu} 81}\end{matrix}$

Comparing Eq. 78-80 to Eq. 19-22 shows the means for the speed andcourse for the two-maneuver case are the same as those for theone-maneuver case, but the variances for the two-maneuver case aresmaller by two thirds. (Recall that τ=T for the one-maneuver case.)

The Man and Variance of the State

Determining the mean and variance of the state when the target undergoestwo (kinematically-constrained) maneuvers is difficult to do inclosed-form. However, these parameters can be approximated using MonteCarlos techniques. For the analysis, 50,000,000 samples were used. Forall the cases, the target was initially at the origin moving due Northat 20 m/s. The update time interval was T=10 s and the maximumacceleration was A=2 m/s². Thus, (x₀, y₀)=(0,0), s₀=20, and θ₀=0. Thetable in FIG. 4 compares the predicted state and covariance from theCGBF against the Kalman filter. Since the Kalman filter does not useturn rates, its predicted state and covariance do not change as themaximum turn rate is varied. However, as the table shows, these valuesvary dramatically depending on the maximum turn rate assumed for thetarget. For example, when the maximum turn rate is small, the CGBFpredicted position is nearly identical to Kalman. But as the maximumturn rate increases to 15 deg/s, the y position diminishes from 200 m to˜151 m, a nearly 25% reduction. The y component of the velocity reducesfrom 20 m/s to ˜8 m/s, which is more than 50% reduction. The varianceterms show an interesting difference as well. The Kalman filter assumesthe distribution is circular; no initial direction of the target isconsidered. As a result, it grossly overestimates the uncertaintyperpendicular to the initial course of the target when it has a smallturn rate. The target needs a turn rate over 10 deg/s for the magnitudeof the uncertainty to be as large as the Kalman filter assumes. Finally,the covariance terms are very different. The Kalman filter onlyspecifies two covariance terms: x position with the x velocity and yposition with the y velocity. The other four covariance terms areassumed to be zero. As shown in the table, this is a reasonableassumption. The last four covariance terms are all near, or possiblyequal to zero regardless of the maximum rate. However, the table showsthat for the two covariance terms that the Kalman filter computes, theyare quite different from the CGBF. The Kalman filters uses ˜667 m²/s forboth x and y. The CGBF estimates the x{dot over (x)} covariance to be assmall as ˜0.2 m²/s for a small maximum turn rate target to as large as˜1565 m²/s for a high maximum turn rate target. Interestingly, the y{dotover (y)} covariance does not vary much relative to the maximum turnrate. For most turn rates, it stays around 400 m²/s, which is about 35%small than the corresponding Kalman value.

The state estimates are directly influenced by the predicted covariance.Sine the predicted covariance from the CGBF is tighter and more alignedto the real target state, it suggests that the CGBF should yield tighterand more accurate state estimates than the Kalman filter.

Motion Update

Performing a motion predict is particularly computational for a GBF.Each cell in the grid must be propagated to all the other cells that thetarget could have transitioned to during the time to the nextmeasurement update. These propagations are found by generating a largenumber of Monte Carlo samples for each cell in the grid. Since the gridstend to be large, an enormous number of Monte Carlo samples are requiredfor each cell to form a good approximation of the cell's transitiondistribution.

To form the target starting position for each cell, there are two commonapproaches: 1) use the center position of the cell, or 2) randomlysample within the cell positional limits. Randomly sampling within thecell limits is usually a better approach because the true targetposition could be near an edge of the cell. If the cells are largecompared to the distance the target moves in an update, then the errorsintroduced by only using the cell center could be significant. However,random sampling over the cell may require even more Monte Carlo sampleper cell to get a good distribution over the entire cell.

With this invention, an alternative motion update procedure has beendeveloped that has two key refinements over these two approaches. Thefirst refinement is how the mass in the cells is transitioned. CurrentGBF approaches move the mass in a cell treating the cell state as a“point”, similar to a particle in a particle filter. In the CGBF,instead of moving a point within the cell, the entire cell region ismoved per Monte Carlo sample. This is illustrated in FIG. 5. By movingthe entire cell for each sample, there is more assurance of obtaining agood distribution over each origination cell, without the need torandomly sample within the cell limits.

Note from FIG. 5 that a transitioned cell from the origination grid willgenerally land across boundaries of cells within the destination grid.The probability in the transitioned origination cell is apportioned outto the destination cells based on the fractional coverage it receives(from the transitioned cell). The function CellOverlapFraction in FIG. 2determines this fractional coverage. The cells do not need to be thesame size from the origination grid to the destination grid. Eachdestination cell only gets the percentage of mass based on the portionof overlap. The second refinement exploits the motion symmetry to obtaintwo Monte Carlo samples from each sample made. This final refinement isdetailed in the next section.

Exploiting the Motion Symmetry

As was pointed out earlier, the transition distribution is symmetricwith respect to the target's initial course. As a result, this symmetrycan be exploited in the generation of the random samples. Each time thetarget is randomly moved throughout the time interval, the mirroredstate can also be used as a sample. Thus, two random samples aregenerated for each random walk which results in 2 m samples beinggenerated for the computational cost of generating m samples.

Let (x₀, y₀) be the (estimated) initial state of the target at thebeginning of the update time interval and (x_(k), y_(k)) be the finalpredicted state at the end of the time interval (by undergoing kmaneuvers). Similarly, let s₀ and θ₀ be the initial speed and course ands_(k) and θ_(k) be the predicted speed and course at the end of theinterval. The mirrored state is defined as the resulting target state ifit had simply reversed all its turns such that left turns become rightturns, and vice versa. This mirrored state can be used as a second MonteCarlo sample for each sample actually computed. The mirrored state isstraightforward to determine directly from the final state of the targetafter it completed all its maneuvers for the time interval.

Referring to FIG. 6, the resulting target state is (x_(k), y_(k)) at theend of time interval. The mirrored target state is (x′_(k), y′_(k)) withspeed s′_(k) and course θ′_(k). Let β be the angle between the initialcourse of the target and the angle to the predicted state, (x_(k),y_(k)). Let α=θ₀+β. Therefore, the mirrored position can be found as:

$\begin{matrix}{{\Delta \; x} = {x_{k} - x_{0}}} & {{Eq}.\mspace{14mu} 82} \\{{\Delta \; y} = {y_{k} - y_{0}}} & {{Eq}.\mspace{14mu} 83} \\{d = \sqrt{{\Delta \; x^{2}} + {\Delta \; y^{2}}}} & {{Eq}.\mspace{14mu} 84} \\{{\sin \; \alpha} = \frac{\Delta \; x}{d}} & {{Eq}.\mspace{14mu} 85} \\{{\cos \; \alpha} = \frac{\Delta \; y}{d}} & {{Eq}.\mspace{14mu} 86} \\{x_{k}^{\prime} = {x_{0} + {d\; {\sin \left( {\theta_{0} - \beta} \right)}}}} & {{Eq}.\mspace{14mu} 87} \\{y_{k}^{\prime} = {y_{0} + {d\; {\cos \left( {\theta_{0} - \beta} \right)}}}} & {{Eq}.\mspace{14mu} 88}\end{matrix}$

But θ₀'2β=θ₀−(θ₀+β)+θ₀=2θ₀−α. So,

$\begin{matrix}\begin{matrix}{x_{k}^{\prime} = {x_{0} + {d\; {\sin \left( {{2\theta_{0}} - \alpha} \right)}}}} \\{= {x_{0} + {d\left\lbrack {{\sin \; 2\theta_{0}\cos \; \alpha} - {\cos \; 2\theta_{0}\sin \; \alpha}} \right\rbrack}}} \\{= {x_{0} + {d\left\lbrack {{\sin \; 2\theta_{0}\frac{\Delta \; y}{d}} - {\cos \; 2\theta_{0}\frac{\Delta \; x}{d}}} \right\rbrack}}} \\{= {x_{0} + {\Delta \; y\; \sin \; 2\theta_{0}} - {\Delta \; x\; \cos \; 2\theta_{0}}}}\end{matrix} & {{Eq}.\mspace{14mu} 89} \\\begin{matrix}{y_{k}^{\prime} = {y_{0} + {d\; {\cos \left( {{2\theta_{0}} - \alpha} \right)}}}} \\{= {y_{0} + {d\left\lbrack {{\cos \; 2\theta_{0}\cos \; \alpha} - {\sin \; 2\theta_{0}\sin \; \alpha}} \right\rbrack}}} \\{= {y_{0} + {d\left\lbrack {{\cos \; 2\theta_{0}\frac{\Delta \; y}{d}} - {\sin \; 2\theta_{0}\frac{\Delta \; x}{d}}} \right\rbrack}}} \\{= {y_{0} + {\Delta \; y\; \cos \; 2\theta_{0}} - {\Delta \; x\; \sin \; 2\theta_{0}}}}\end{matrix} & {{Eq}.\mspace{14mu} 90}\end{matrix}$

Computing the mirrored speed and course is starightforward as well. Bydefinition, the mirrored speed does not change; only the target'scourse. Let Δθ=θ_(k)'1θ₀ be the course difference from the initialcourse to the final course. Thus, the final course can be thought of asθ_(k)=θ₀+Δθ. The mirrored course must be θ_(k)=θ₀ 31 Δθ. Therefore themirrored state is:

x ′ _(k) =x ₀+(y _(k) −y ₀)sin 2θ₀−(x _(k) −x ₀)cos 2θ₀   Eq. 91

y′ _(k) =y ₀+(y _(k) −y ₀)cos 7θ₀+(x _(k) −x ₀)sin 2θ₀   Eq. 93

s′_(k)=s_(k)   Eq. 93

θ′_(k)=θ₀−Δθ=θ₀−(θ_(k)−θ₀)=2θ₀−θ_(k)   Eq. -b 94

Using Uniform RVs to Approximate Normal RVs

Although many computer systems provide the ability to generate normalrandom variables, they typically require many more computations than thegeneration of uniform random variables. When performing the Monte Carlosampling, the random variables must be drawn from a normal distribution(see FIG. 2). Fortunately, the CGBF sampling process does not need“true” normal random variables; they only need to be approximatelynormal. A simple approximation was developed for the CGBF to exploitthis idea. The normal-like random variables are generated using thefollowing procedure. Suppose a total of m random samples are neededwhere each random sample is drawn from a normal distribution with aspecified standard deviation. The number of random samples is dividedinto thirds. Let k=m/3. Then:

-   -   generate k uniform random samples that are within one sigma,    -   generate k uniform random samples that are within two signma,        and finally,    -   generate k uniform random samples that are within three sigma.

This method is called the thirds procedure. For a “true” (1D) normalrandom variable, ˜68.2% of the samples would be within one sigma, ˜95.4%are within two sigma, and ˜99.7% are within three sigma[12]. Using theapproximating procedure just described,

${k + {\frac{1}{2}k} + {\frac{1}{3}k}} = {{\frac{k}{6}\left( {6 + 3 + 2} \right)} = {{\frac{11}{6}\left( \frac{m}{3} \right)} = {\left( {61.1\%} \right)m}}}$

samples will be within one sigma,

$k = {{k + {\frac{2}{3}k}} = {{\frac{k}{3}\left( {3 + 3 + 2} \right)} = {{\frac{8}{3}\left( \frac{m}{3} \right)} = {\left( {88.8\%} \right)m}}}}$

samples will be within two sigma, and of course, 100% of the samples wilbe within three sigma. These percentages are close enough to get theproper sampling. Thus, using this thirds procedure with uniform randomvariables, the needed normal random variables are efficiently obtained.

Although the invention has been described in detail with particularreference to these preferred embodiments, other embodiments can achievethe same results. Variations and modifications of the present inventionwill be obvious to those skilled in the art and it is the indent of thisapplication to cover, in the appended claims, all such modifications andequivalents. The entire disclosure and all references, applications,patents and publications cited above are hereby incorporated byreference.

What is claimed is:
 1. A computer program product embodied in anon-transitory computer-readable medium that provides state estimatesfor a maneuvering target using a multi-dimensional grid-based filter,comprising: software instructions for enabling a computer to performpredetermined operations; and a machine-readable medium bearing thesoftware instructions, the predetermined operations comprising: assuminggiven target kinematic constraints, forming a first grid of cells, of afinite size and scale, wherein each cell corresponds to a particularrectangular region in a x and y coordinate system, capturing possibletarget states using a first measurement, wherein each state is aposition, course and speed of the target; forming a second grid of cellsof finite size and scale, wherein each cell corresponds to a particularrectangular region in an x and y coordinate to over a duration betweenthe first measurement and a second measurement; forming a predictedtarget state distribution in the second grid by sampling the possibletarget states represented in the first grid, wherein the predictedtarget state distribution is represented in the cells of the secondgrid, wherein each cell includes: a probability that the targettransitioned to the cell; a mean speed that the target would have if ittransitioned to the cell; a speed variance that the target would have ifit transitioned to the cell, wherein the speeds within each cell areassumed normally distributed having the mean speed and speed variance; amean course that the target would have if it transitioned to thecorresponding rectangular region; and a course variance that the targetwould have if it transitioned to the cell, wherein courses within eachcell are assumed normally distributed having the mean course and coursevariance; and forming a discretized probability distribution over allthe cells in the second grid; and substituting the second grid as firstgrid for a subsequent measurement.
 2. The computer program product ofclaim 1, wherein the finite size and scale of each grid is based on ameasurement error distribution associated with each grid, wherein thefinite size of the grid is 20×20 cells or less.
 3. The computer programproduct of claim 1, wherein the sampling simulates a different randommaneuver of the target over a specified duration, and symmetry of thepredicted state distribution is used to double the effect of thesampling.
 4. The computer program product of claim 1, wherein cell-sizedparticles are transitioned from the first grid to the second grid duringthe sampling and the cell probability of the transitioned cell from thefirst grid is apportioned to the cells it overlaps in the second grid.5. A method for providing state estimates for a maneuvering target usinga multi-dimensional grid-based filter, comprising: assuming given targetkinematic constraints, forming a first grid of cells, of a finite sizeand scale, wherein each cell corresponds to a particular rectangularregion in an x and y coordinate system, capturing possible target statesusing a first measurement, wherein each state is a position, course andspeed of the target; forming a second grid of cells of finite size andscale, wherein each cell corresponds to a particular rectangular regionin an x and y coordinate system, comprising possible positions thetarget could have transitioned to over a duration between the firstmeasurement and a second measurement; forming a predicted target statedistribution in the second grid by sampling the possible target statesrepresented in the first grid, wherein the predicted target statedistribution is represented in the cells of the second grid, whereineach cell includes: a probability that the target transitioned to thecell; a mean speed that the target would have if it transitioned to thecell; a speed variance that the target would have if it transitioned tothe cell, wherein the speeds within each cell are assumed normallydistributed having the mean speed and speed variance; a mean course thatthe target would have if it transitioned to the correspondingrectangular region; and a course variance that the target would have ifit transitioned to the cell, wherein courses within each cell areassumed normally distributed having the mean course and course variance;and forming a discretized probability distribution over all the cells inthe second grid; and substituting the second grid as a first grid for asubsequent measurement.
 6. The method of claim 5, wherein the finitesize and scale of each grid is based on a measurement error distributingassociated with each grid, wherein the finite size of the grid is 20×20cells or less.
 7. The method of claim 5, wherein the sampling simulatesa different random maneuver of the target over a specified duration, andsymmetry of the predicted state distribution is used to double theeffect of the sampling.
 8. The method of claim 5, wherein cell-sizedparticles are transitioned from the first grid to the second grid duringthe sampling and the cell probability of the transitioned cell from thefirst grid is apportioned to the cells it overlaps in the second grid.